% % File		:	verity/ppg_0_1_display
% % Description	:   for Heart Rate Algorithm
% % Fuction		:	PPG signal display
% % Author		:	@hiyangdong
% % Version		:	V0.1
% % Time		:	12 Jan. 2016

% 0.1.1 read data and orignal calculate
close all
clear all
orignal_txt = textread('PAH0119_1022.txt','%s');    %#ok<DTXTRD> %read txt file data 32bit with 4byte and Low byte first
orignal_dec = hex2dec(orignal_txt);                 %hex 2 dec every byte
orignal_num = length(orignal_dec);                  %length of the orignal data
% orignal_num = 1000;                  %length of the orignal data
ppg_data_byte = 20;
ppg_num  = orignal_num/ppg_data_byte;               %length of the PPG data
ppg_data = zeros(1,ppg_num);                        %generate an 1-by-length empty array

% % 0.1.2 convert orignal hexadecimal data to ppg 32bit array
% j=1;
% for i=1:4:orignal_num
%     ppg_data(j)=(orignal_dec(i)+(orignal_dec(i+1)*256)+(orignal_dec(i+2)*256*256)+(orignal_dec(i+3)*256*256*256));
%     j=j+1;
% end

% % 0.9.1 convert orignal hexadecimal data
% % data format:
% % 1 byte address
% % 2 byte PPG orignal data
% % 6 byte delta ppg data relative to it's fore-data --9 byte
% % 1 byte Gyro X
% % 1 byte Gyro Y
% % 1 byte Gyro Z
% % 1 byte Acce X
% % 1 byte Acce Y
% % 1 byte Acce Z
% % 1 byte Magn X
% % 1 byte Magn Y
% % 1 byte Magn Z --9 byte
% % 1 byte Barometer Z
% % 1 byte Temprature --2 byte
%
% j=1;
% for i=1:ppg_data_byte:orignal_num
%     ppg_data(j)=orignal_dec(i+1)*256 + orignal_dec(i+2);
%     j=j+1;
%     ppg_data(j)=orignal_dec(i+3)*256 + orignal_dec(i+4);
%     j=j+1;
%     ppg_data(j)=orignal_dec(i+5)*256 + orignal_dec(i+6);
%     j=j+1;
%     ppg_data(j)=orignal_dec(i+7)*256 + orignal_dec(i+8);
%     j=j+1;
% end
%
% % 0.9.2 convert orignal hexadecimal data to ppg 16bit array
% j=1;
% m = 9;
% for i=1:ppg_data_byte:orignal_num
%     gyro_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     gyro_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     gyro_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     acce_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     acce_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     acce_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     magn_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     magn_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     magn_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     baro_data(j) = orignal_dec(i+m);
%     j=j+1;
% end
%
% j=1;
% m = m+1;
% for i=1:ppg_data_byte:orignal_num
%     temp_data(j) = single((orignal_dec(i+m)+100)/10.0);
%     j=j+1;
% end

% % 1.1.1 convert orignal hexadecimal data
% data format first frame:
% 1 byte address
% 2 bytes PPG orignal data
% 2 bytes Gyro X
% 2 bytes Gyro Y
% 2 bytes Gyro Z
% 2 bytes Acce X
% 2 bytes Acce Y
% 2 bytes Acce Z

% data format second frame:
% 1 byte address
% 2 bytes PPG orignal data
% 2 bytes Gyro X
% 2 bytes Gyro Y
% 2 bytes Gyro Z
% 2 bytes Magn X
% 2 bytes Magn Y
% 2 bytes Magn Z
% 2 bytes Barometer Z
% 2 bytes Temprature

j=1;
for i=1:ppg_data_byte:orignal_num
    ppg_data(j)=orignal_dec(i+1)*256 + orignal_dec(i+2);
    j=j+1;
end

% 1.1.2 convert orignal hexadecimal data to ppg 16bit array
j=1;
m = 3;
for i=1:ppg_data_byte:orignal_num
    temp_data = orignal_dec(i+m)*256 + orignal_dec(i+m+1);
    gyro_x_data(j) = typecast (uint16(temp_data),'int16');
    j=j+1;
end

j=1;
m = 5;
for i=1:ppg_data_byte:orignal_num
    temp_data = orignal_dec(i+m)*256 + orignal_dec(i+m+1);
    gyro_y_data(j) = typecast (uint16(temp_data),'int16');
    j=j+1;
end

j=1;
m = 7;
for i=1:ppg_data_byte:orignal_num
    temp_data = orignal_dec(i+m)*256 + orignal_dec(i+m+1);
    gyro_z_data(j) = typecast (uint16(temp_data),'int16');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    acce_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_x_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_y_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    magn_z_data(j) = typecast (uint8(orignal_dec(i+m)),'int8');
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    baro_data(j) = orignal_dec(i+m);
    j=j+1;
end

j=1;
m = m+1;
for i=1:ppg_data_byte:orignal_num
    temp_data(j) = single((orignal_dec(i+m)+100)/10.0);
    j=j+1;
end

% % 0.4.2 ppg data max and min analysis
% ppg_data_max = max(ppg_data)
% ppg_data_min = min(ppg_data)
% ppg_data_minus = ppg_data_max - ppg_data_min
% for i=1:ppg_num-1
%     ppg_data_inc(j)=ppg_data(i+1) - ppg_data(i);
%     j=j+1;
% end
% ppg_data_inc_max = max(ppg_data_inc)
% ppg_data_inc_min = min(ppg_data_inc)
% ppg_data_inc_minus = ppg_data_inc_max - ppg_data_inc_min
% num232 = power(2,32)
% num224 = power(2,24)
% num216 = power(2,16)
% ppg_data = ppg_data(1:length(ppg_data)/4);
% ppg_data(length(ppg_data)/4:2*length(ppg_data)/4) = mean(ppg_data(length(ppg_data)/4:2*length(ppg_data)/4));
% 0.1.3 time domain analysis
ppg_num = length(ppg_data);
sampletime = 0.01;          %sample time that unit is second
t = 0:sampletime:(ppg_num-1) * sampletime;      %time point

% 0.1.4 time domain plot
figure(1);
clf();
subplot(2,2,1);
plot(t,ppg_data,'k');
hold on
title('Time domain orignal PPG data');

% 0.1.5 frequency domain basic analysis
samplefreq=1/sampletime;    %sample frequency that unit is Hz
f=samplefreq*(0:(ppg_num-1))/ppg_num;           %frequency point

% 0.1.6 frequency domain offset calculate analysis
% sampledata_ave=0;
% for i=1:ppg_num
%     sampledata_ave = sampledata_ave+ppg_data(i);
% end
% sampledata_ave=sampledata_ave/ppg_num;
% sampledata_offset=ppg_data-sampledata_ave;

% 0.1.6.1 modify the average algorithm
sampledata_offset=ppg_data-mean(ppg_data);

% 0.1.6.2 add the plot of offset
% sampledata_offset2=sampledata_offset/30+127;
% plot(t,sampledata_offset2,'k');

% 0.1.7 frequency domain FFT analysis
fft_sampledata=fft(sampledata_offset);
fft_sampledata_abs=abs(fft_sampledata);

% 0.1.8 frequency domain plot
N=ppg_num;
subplot(2,2,3)
plot(f(1:N/2),fft_sampledata_abs(1:N/2),'k');
title('Frequency domain orignal PPG data ');

% 0.3.1 lowpass filter design
fpass = 5
fstop = samplefreq;
f_num = fix(fpass*length(t)/samplefreq);
fft_sampledata(f_num:(length(t)-f_num+1)) = 0;

% 0.3.2 lowpass filtered data
ppg_fliterdata = real(ifft(fft_sampledata)) + mean(ppg_data);

% 0.3.3 lowpass filtered data display
subplot(2,2,2);
plot(t,ppg_fliterdata);
hold on
title('Time domain filtered PPG data');

% 0.3.4 filtered data frequency domain FFT analysis
fft_fliterdata=fft(ppg_fliterdata);
fft_fliterdata_abs=abs(fft_fliterdata);

% % 0.3.5 filtered data frequency domain plot
% subplot(2,1,2);
% plot(f(1:N/2),fft_fliterdata_abs(1:N/2),'b');
% title('Frequency domain filtered PPG data');

% 0.3.5.1 filtered data frequency domain plot
subplot(2,2,4);
fft_fliterdata_abs = abs(fft_sampledata);
plot(f(1:N/2),abs(fft_sampledata(1:N/2)),'b');
title('Frequency domain filtered PPG data');

% 0.2.1 peaks and throughs detection
ds=diff(ppg_fliterdata);
filter=find(ds(2:end)==0)+1;%%find zeros
ds(filter)=ds(filter-1);%%replace zeros
ds=sign(ds);
ds=diff(ds);
throughs=find(ds>0); %when 0 change to 2 the point disappear
peaks=find(ds<0);

% 0.2.2 peaks and throughs addtional plot
subplot(2,2,2);
plot(t(peaks+1),ppg_fliterdata(peaks+1),'or');
plot(t(throughs+1),ppg_fliterdata(throughs+1),'b*');

% % 0.5.1 heart beat calculate
% ppg_time = max(t)  % t is no problem which I check with my watch
% ppg_peak = length(peaks)
% heart_rate = ppg_peak/ppg_time*60

% 0.5.1.1  heart beat calculate
ppg_peak = 0;
ppg_time = max(t)
peaks_num = length(peaks)
for k = 1:(peaks_num-1)
    if (t(peaks(k+1)) - t(peaks(k)))>0.3
        ppg_peak = ppg_peak+1;
    end
end
ppg_peak;
heart_rate = ppg_peak/ppg_time*60

% % 0.8.3  heart beat calculate
% figure(2);
% subplot(2,3,5);
% plot(t(1:N/4),gyro_x_data,'k');
% % set(gca,'YLim',[min(gyro_x_data)-10,max(gyro_x_data)+10]);
% subplot(2,3,6);
% plot(t(1:N/4),gyro_z_data,'k');
% set(gca,'YLim',[min(gyro_z_data)-10,max(gyro_z_data)+10]);

% % 1.0.0  heart beat calculate
figure(2);
plot1x = gyro_x_data;
plot1y = gyro_y_data;
plot1z = gyro_z_data;
plot2x = acce_x_data;
plot2y = acce_y_data;
plot2z = acce_z_data;

subplot(2,3,1);
plot(t,plot1x,'k');
scale = (max(plot1x) - min(plot1x))/5;
set(gca,'YLim',[min(plot1x)-scale,max(plot1x)+scale]);

subplot(2,3,2);
plot(t,plot1y,'k');
scale = (max(plot1y) - min(plot1y))/5;
set(gca,'YLim',[min(plot1y)-scale,max(plot1y)+scale]);

subplot(2,3,3);
plot(t,plot1z,'k');
scale = (max(plot1z) - min(plot1z))/5 + 1;
set(gca,'YLim',[min(plot1z)-scale,max(plot1z)+scale]);

subplot(2,3,4);
plot(t(1:N/4),plot2x,'k');
scale = (max(plot2x) - min(plot2x))/5;
set(gca,'YLim',[min(plot2x)-scale,max(plot2x)+scale]);

subplot(2,3,5);
plot(t(1:N/4),plot2y,'k');
scale = (max(plot2y) - min(plot2y))/5;
set(gca,'YLim',[min(plot2y)-scale,max(plot2y)+scale]);

subplot(2,3,6);
plot(t(1:N/4),plot2z,'k');
scale = (max(plot2z) - min(plot2z))/5;
set(gca,'YLim',[min(plot2z)-scale,max(plot2z)+scale]);

% % 1.0.1  heart beat calculate
figure(3);
plot3x = magn_x_data;
plot3y = magn_y_data;
plot3z = magn_z_data;

subplot(1,3,1);
plot(t(1:N/4),plot3x,'k');
scale = (max(plot3x) - min(plot3x))/5;
set(gca,'YLim',[min(plot3x)-scale,max(plot3x)+scale]);

subplot(1,3,2);
plot(t(1:N/4),plot3y,'k');
scale = (max(plot3y) - min(plot3y))/5;
set(gca,'YLim',[min(plot3y)-scale,max(plot3y)+scale]);

subplot(1,3,3);
plot(t(1:N/4),plot3z,'k');
scale = (max(plot3z) - min(plot3z))/5 + 1;
set(gca,'YLim',[min(plot3z)-scale,max(plot3z)+scale]);

% 1.0.2
figure(3);
plot4_1 = baro_data;
plot4_2 = temp_data;

subplot(1,2,1);
plot(t(1:N/4),plot4_1,'k');
scale = (max(plot4_1) - min(plot4_1))/5;
set(gca,'YLim',[min(plot4_1)-scale,max(plot4_1)+scale]);

subplot(1,2,2);
plot(t(1:N/4),plot4_2,'k');
scale = (max(plot4_2) - min(plot4_2))/5;
set(gca,'YLim',[min(plot4_2)-scale,max(plot4_2)+scale]);
